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Starting from the second equilibrium equation in the BBGKY hierarchy under the 
Kirkwood superposition closure, we implement a new method for studying the asymp- 
totic decay of correlations in the hard disk fluid in the high density regime. From our 
analysis and complementary numerical studies, we find that exponentially damped 
oscillations can occur only up to a packing fraction rj* ~ 0.718, a value which is in 
substantial agreement with the packing fraction, rj ~ 0.723, believed to characterize 
the transition from the ordered solid phase to a dense fluid phase, as inferred from 
Mak's Monte Carlo simulations [Phys. Rev. E 73, 065104 (2006)]. We next show 
that the same method of analysis predicts that exponential damping of oscillations 
in the hard sphere fluid becomes impossible when A = 4n7Rx 3 [l + H(l)] > 34.81, 
where H(l) is the contact value of the correlation function, n is the number density 
and a is the sphere diameter, in exact agreement with the condition, A > 34.8, first 
reported in a numerical study of the Kirkwood equation by Kirkwood et al. [J. Chem. 
Phys. 18, 1040 (1950)]. Finally, we show that our method confirms the absence of 
any structural transition in hard rods for the entire range of densities below close 
packing. 
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I. INTRODUCTION 



The second equilibrium equation in the BBGKY hierarchy establishes an exact relation 
between the pair and triplet number density. Invoking the Kirkwood superposition approx- 
imation yields a nonlinear integral equation for the pair correlation function^. Interest in 
studying the analytic and numerical properties of the resulting Yvon-Born-Green and/or 
Kirkwood equation began with Kirkwood and coworkers^-, and continues to the present 
day 7 . Of particular interest is whether the closed equation provides an essentially correct 
description of the fluid phase, and whether a (possible) change in the analytic character 
of the solutions signals a change from the fluid phase to a solid phase. It is to the latter 
question that the methods of the present contribution are directed. 

One approach to explore analytically the possibility of a phase transition from the fluid 
phase to the solid phase is to mobilize the theory of nonlinear integral equations, focusing 
on theorems which establish the necessary and sufficient conditions for the existence and 
uniqueness of solutions, and bifurcation points^ - —. An alternative approach is to introduce 
a moment expansion by means of which the YBG equation can be cast into a nonlinear 
differential equation which may be used to analyze long-range correlations^ 1 ^. The present 
contribution is centered on a new method of studying the asymptotic decay of correlations, 
first introduced in Ref. 15 for the hard sphere fluid. 

In this paper, we focus on the hard disk fluid. The method leads to the prediction 
of a structural transition in both the hard sphere and hard disk fluids, and no transition 
in the hard rod system (as must be the case). We shall show that the values of packing 
fractions at which the predicted transitions occur are in agreement with estimates derived 
from numerical solution of the Kirkwood equation^ and recent Monte Carlo simulations^. 

II. THE SECOND EQUILIBRIUM HIERARCHY EQUATION FOR HARD 
DISKS 

We consider a gas of hard disks of diameter a at thermal equilirium with constant number 
density n and temperature T. 

Let n 2 (r 12) denote the number density of pairs of particles situated at distance r\ 2 = 
|ri — r 2 |. Using the fact that n 2 (r 12 ) = for r 12 < a we introduce a dimensionless function 
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2/2(^12) defined by 

n 2 (r 12 ) = n 2 9(r 12 - a)y 2 (r 12 ) (1) 

where 9 is a unit step function. 

We assume that n 2 (r 12 ) — > n 2 when r 12 — > 00. The two-particle dimensionless correlation 
function h 2 (ri 2 ) is then defined by the cluster decomposition y 2 (r\ 2 ) = 1 + h 2 (ri 2 ), so that 
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n 2 



0(ru - a)[l + h 2 (r 12 )] (2) 



h 2 (ri 2 ) is thus supposed to satisfy the asymptotic condition 

lim h 2 (r 12 ) = (3) 

ri2— >oo 

The second equilibrium Yvon-Born- Green (YBG) hierarchy equation establishes an ex- 
act relation between n 2 (ri 2 ) and the reduced three-particle number density n^(ri,r 2 ,r^). 
Introducing the excluded volume factor we write the three-particle density as 

ri 3 (ri, r 2 , r 3 ) = 9(r 12 - a)9(r 13 - a)9(r 23 - a)n 3 y 3 (r 12 , r 13 , r 23 ) (4) 

The function n 3 depends only on the distances = |r^-| = |r^ — r,'|, and is a symmetric 
function of the three variables ri2, ?"i3, r 23 . In the case of hard disks y 2 is related to y 3 through 
the second YBG hierarchy equation (see Appendix A) 

d f 

—2/2(7") = na / d&(r ■ &)9(\r - a&\ - a)y 3 (r, a, |r - aa\) (5) 

where <x and r are unit vectors |<x| = |f| = 1. Putting r ■ & = cos0 we rewrite ([5]) in an 
explicit form 

-j-y 2 (r) = na / d(p cos y 3 (r, a, \/r 2 — 2ra cos + a 2 )9{r — 2a cos 0) (6) 
When writing fl6]) the equality 

#(|r — acr\ — a) = 9{r — 2a cos 0) 

has been used. 

The rigorous relation fl6]), valid for r > a, will be the starting point for subsequent 
considerations. 
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III. THE KIRKWOOD SUPERPOSITION APPROXIMATION 



The Kirkwood superposition approximation consists in replacing in equation the 
three-particle density by the product of two-particle densities corresponding to three different 
pairs of particles 

y 3 (r,a,\r-a&\) ->■ y 2 (r)y 2 ((r)y 2 (\r - aa\) (7) 
Adopting (j7|) leads to a closed equation 

~i~y2(f) = nay2(r)y2(cr) I d<p cos 2/2 (V r2 _ 2r<xcos0 + a 2 )9{r — 2crcos0) (8) 
" r Jo 

It is convenient to rewrite (|SJ) using the dimensionless distance x — r/a. Denoting by Y(x) 

the function 

= y 2 (xa) (9) 
we find that it satisfies the non-linear equation 

— lnF(x) = na 2 Y{\) I d<p cos0 Y( x 2 - 2x cos + 1)6 (x - 2 cos0) (10) 
ax Jo 

valid in the region 2 > 1. Equation ( flOl) represents the closure of the YBG hierarchy 
corresponding to the superposition approximation. 

Our aim is to derive from (fTOl) the equivalent integral equation satisfied by the dimen- 
sionless correlation function 

H(x) = Y(x)-l (11) 
To this end we insert fflTj) into ffTOl finding 

— ln[l + = na 2 [l + if (1)] <^ / d0cos0#(x - 2cos0) (12) 

dx I ./n 



+ J d(f> cos (f)H (a/x 2 — 2a; cos0 + l)0(a; — 2cos0)| 

In order to derive an integral equation for H(x) we integrate ffl2l over the spatial interval 
(x, 00). Using the formulae 



J d(f) 6 (x- 2 cos (f)) cos (f)= -20(2 (|) 2 (13) 

(14) 



y dz6{2-z)yjl- (£)' = 0(2- x) -fyl - (|) 2 + arccos | 
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we get 



where 



and 



ln[l + H(x)} = 2na 2 [l + H{l)][I*{x) + I**{x)\ 



I*(x) = 9(2 -x) 



x 



X 

, + arccos — 
2) 2 



poo p2tt 

I**(x) = — / dz dd> cos 4>6(z - 2cos0)#(v / ^ 2 + 1 - 2^cos< 
2 Jx Jo 



(15) 



(16) 



It turns out that the angular integration in the formula for I**(x) can be explicitly performed 
(the calculation is presented in Appendix B). One eventually finds 



f x+1 fx 2 + s 2 - 1 

I**(x) = — / ds 6(s — 1) sH(s) arccos 

Jx-l 



2xs 



(17) 



In this way we arrive at an integral equation 



ln[l + H(x)} = 2na 2 [l + H(l)) { 6(2 - x) 



x 



x 



\l 1 — I — ) + arccos 

2V \2J 2 



18) 



x+l / x 2 s 2 _ ]_ 

ds 9(s — 1) sH(s) arccos 

x-l 



2xs 



representing the superposition closure for the correlation function H(x). 

Our method of determining H(x) is based on the fact that the solution of Eq. (1181) satis- 
fying the boundary condition lim^oo H(x) = can be obtained numerically by iterations. 
Note that we can rewrite Eq. f|T8|) in the following form 



H(x) = C(H(x)) 



(19) 



where £ is the integral operator given by 



C(H(x)) = -1 + exp 2na 2 [l + H(l)] < 6(2 - x) 



-fv 1 -^ 



x 

arccos — 
2 



x+l 



ds6(s — 1) sH(s) arccos 



x-l 



x 



2 ^s 2 -l 



2xs 



(20) 



The above integral equation for H(x) was solved by a standard Neumann method with 
succesive over-relaxation^. The iterative solutions are then given by 



H n = (1 - a)H n ^ + aC(H n ^) 



(21) 
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FIG. 1. Pair correlation function H(x) for the surface fraction £ = ra7rcr 2 /4 = 0.35 (dashed line) 
and £ = 0.59 (solid line). 



The relaxation parameter a was taken to be 0.25. The iterations were continued until 
successive values of H(x = 0) differed by less than e = 10 -5 , except in the vicinity of 
the threshold surface fraction £* (see Sec. IV Bj) . where the convergence was slow and the 
iterations were discontinued at e = 10~ 2 . 

Examples of the correlation functions obtained in this way are presented in Fig. [IJ As 
it is seen, the decay of H(x) becomes slower as the surface fraction is increased, and a 
pronounced peak structure appears. 



IV. LINEARIZATION OF KIRKWOOD'S EQUATION: THE RING 
APPROXIMATION 

Before continuing the analysis based on equation (1181) let us make a comment on the rela- 
tionship between the superposition approximation and the ring approximation, well known 
from the kinetic theory (see^ and references given therein). 

Originally, the ring approximation was applied to the study of long wavelength hydrody- 
namic phenomena, and was defined by neglecting in the second equation of the dynamical 
BBGKY hierarchy all contributions from the three-particle correlations. The three-particle 
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correlation function h 3 is defined by the cluster decomposition of the density y 3 



2/3(^12, ri 3 ,r 2 3) = 1 + h 2 {r 12 ) + h 2 (r 13 ) + h 2 (r 23 ) + h 3 (r 12 ,r 13 ,r 23 ) 



Neglecting h 3 is thus equivalent to the approximation 



2/3(712, ria, r 23 ) ~ 1 + /i2(yi2) + h 2 (r 13 ) + h 2 (r 23 ) 



(22) 



(23) 



It is important here to note that while rejecting in the cluster decomposition (|22|) three- 
particle correlations we nevertheless retain in the hierarchy equation the full excluded volume 
factor represented by the product of unit step- functions (see definition (j3J). In fact, this 
factor represents the exact lowest order term in the density expansion of the three-particle 
number density, and is of fundamental importance for the correct description of hard disks 
at low densities. 

Comparing fl2"31 with the superposition formula 



2/3(^12, ri 3 ,r 23 ) ~ [1 + h 2 (r 12 )\ [1 + h 2 (r 13 )][l + h 2 (r 23 )) 



(24) 



we see that the ring approximation corresponds exactly to the linearization of the Kirkwood 
theory in h 2 . 

The linearized form of equation (j!2p reads 



d 

dx 



H{x) = no 2 I -20(2 - x)^l- (|) 2 [1 + H(l) + H{x)\ 



(25) 



+ j dip cos(p9(x — 2 cos (j))H(\/ ' x 2 + 1 — 2xcos0) | 

In order to derive an integral equation for H(x) we proceed as before by integrating (125 
over the spatial interval (x, 00). The result reads 

H(x) = 2na 2 [J*{x) + J**(x)} 

where 
J*(x) = 6(2 -x) 



x 



— A/1— — + arccos — 
2V \2J 2 



[1 + H(1)]+J dzyjl- ( 



(26) 



and 



J**( x ) = — / dz 



COS 0O[Z 



2 COS (p)H(yJz 2 + 1 - 2z COS . 



r*(x) 



(see equation ([TB]) ). Using again the calculation presented in Appendix B we eventually find 



f x+1 fx 2 + s 2 -l' 
J**( x ) = - I ds 6(s - 1) sH(s) arccos ) (27) 

Jx-l V ^ XS 



In this way we arrive at an integral equation 
H(x) = 2na 2 6(2-x) 



x I (x\ x 
1 — — + arccos 



2 V \2J 2 



[l + H(l)] + j\z^l-Q 2 H(z)] 

f x+1 fx 2 + s 2 -l\ 
2na 2 l ds 9{s — 1) sH(s) arccos I ] 



(28) 



representing the ring approximation. Equation f[28|) can again be solved numerically by 
iterations, and is expected to yield relevant results in the low density regime. 

The comparison of the results obtained with the full Kirkwood approximation and its 
linearization is presented in Fig. |2j Defining the surface fraction £ as 

9 

£ = ™ U - (29) 

we plot here the compressibility factor Z(£), given by 

Z(0 = = 1 + \no 2 (\ + H(l)) = 1 + 2£(1 + H(l)), (30) 

where the contact value H(l) is also a function of £. For comparison, we include here Z(£) 
dependence as predicted by the scaled particle theory (SPT) 19 

ZspAO = (31) 

As it is seen, the iteration results agree with the SPT predictions up to approximately 
£ = 0.4. For larger packing fractions, the Kirkwood approximation tends to underestimate 
the compressibility factor, whereas ring approximation overestimates it. 



V. ASYMPTOTIC DECAY OF CORRELATIONS: PREDICTING A 
STRUCTURAL TRANSITION 

A. Breakdown of the method used for attractive interactions 

The fundamental information concerning the internal structure of the system is contained 
in the spatial dependence of correlations. In particular, the law governing the asymptotic 
vanishing of correlations is of primary importance. 
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FIG. 2. Compressibility factor Z versus the surface fraction £ for the full Kirkwood approximation 
(circles), the ring approximation (squares), and the scaled particle theory (Eq. (f3T|) . dashed line) 

In order to determine the behavior of H(x) at large distances, it seems natural to follow 
the moment analysis presented in Refs. 13 and 14. The calculation would proceed as follows. 
For x > 2, the Kirkwod equation ( fT2l) can be conveniently written as 

^-\n[l + H(x)]=na 2 [l + H(l)} ! da{± ■ ar)H{\x - <x|) (32) 

When x 3> 1, the power series expansion of -ff(|x — &\) around the point |x| = x yields 
nonzero contributions only from terms involving odd powers of cos0 = (x ■ &). The calcu- 
lation up to the third derivative of H yields the expansion 

[ln[l + H(x)]' = na 2 [1 + H(l)} / cos0{- cos <f)H'{x) (33) 

Jo 

_^ [COS(/) _( COS 0)3] -I(COS0) 3 F'"(X) + ...} 

where ' denotes the derivative with respect to x. 

Using the boundary condition lim x ._ )>00 H(x) = together with the equalities 

d(j)(cOS (j)) 2 = 71, / G?0(COS0) 4 = —71, 

Jo 4 
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and adopting for large distances the asymptotic formula 



\n[l + H(x)] ~H{x) 



we find a linear differential equation of the form 



H"{x) + -H'(x) + a 2 H = 



(34) 



x 



where 



a 2 = 8(1 + 



1 



)>0 



*[1+H(1)) 



But (134]) is the equation for the Bessel function J (ax). The solution of (|34|) is thus an 
oscillating function whose amplitude decays as 1/ y/x . 

Unfortunately, the above result is in disagreement with numerical predictions of expo- 
nentially damped oscillations. Use of the moment expansion as developed in Ref. 13 leads 
to erroneous results. Clearly one has to take into account the whole infinite series in the 
expansion of H(\x — <x|) to get reliable predictions. We have thus to give up this kind of 
expansion, and look for a different approach. 

In broad outline, the failure of the moment expansion developed in Refs. 13 and 14 to 
describe the structural transition in the hard disk fluid can be traced to the nature of the 
governing intermolecular potential. The studies^ 1 ^ deal with the description of correlations 
in the vicinity of the liquid- vapor critical point, where attractive forces play a dominant role. 
Using the moment expansion, one recovers at large distances the classical Ornstein-Zernike 
formula. The present study deals with the fusion transition, where short-range, repulsive 
forces play the critical role. The new method introduced here (see Sec. IV Bl) effectively 
accounts for the difference in the potential governing these two transitions and, anticipating 
our later results, leads to an analytic prediction of the packing fraction at which a structural 
transition occurs in the hard disk (and hard sphere) fluid. 

B. Prediction of a structural transition 

Let us consider again the region x > 2 where the equation (1321) holds. As lim^oo H(x) = 
0, we can replace in (132]) the function ln[l + H(x)} by H(x), and consider the equation 




(35) 
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where 



A = 2ixno 2 [1 + H(l)) 



We then use the expansion 



sm 



2 



\/ x 1 — 2x COS 0+1 = X — COS (j) + 



+ ... 



(36) 



to arrive at the equation 



d 

dx 




(37) 



valid for x ^> 1. 

In order to determine the large x behavior of correlations we have thus to analyze the 
solution of (|37|) . We notice that all derivatives of H(x) satisfy the same equation. It is thus 
natural to consider H( linear combination of exponential modes exp(nx), where k is 

a complex number. The function exp(/cc) satisfies ( 13T1) provided k solves the equation 



where is a modified Bessel function. The physically acceptable solutions k = a + ib are 
those with negative real part a < which assures exponential damping of oscillations. The 
first root (with the smallest absolute value of a), corresponding to the slowest decay of the 
correlation function is shown in Fig. [3J The values of k(£) predicted with the use of Eq. (13"8"j) 
are in good agreement with the decay of the amplitude of H(x) determined by the iterative 
solution of integral equation (|T9|) . Plotting the absolute value of H(x) on a logarithmic plot, 
and fitting it by the single mode ae K ' x , we obtain values of k' which are slightly below those 
obtained by solving Eq. (|38|) . For example, for £ = 0.3 we get k' ~ —2.1 (c/. Fig.H]), whereas 
the corresponding value of k for that surface fraction is k ~ —2.02. Similarly, for £ = 0.55 
we get respectively k' ~ —0.8 and k ~ —0.7. The fact that the values of k' remain slightly 
below those of k can be understood by noting that k corresponds to the slowest decaying 
mode, whereas in the numerical data on H(x) we also see nonzero contributions from other, 
faster decaying modes. 

An interesting feature of the «(£) dependence presented in Fig. [3] is the fact that the 
root becomes purely imaginary at A = A* 15.1. This point can be made more precise 
by analyzing when the equation (I3"8l acquires a purely imaginary solution k = ib. As 




Ah(K) 



(38) 
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FIG. 3. The real part (solid line) and imaginary part (dashed line) of the root of Eq. (|38p corre- 
sponding to the slowest decaying mode. 



I\{ib) = iJ\(b), (1381) implies the condition 

Mb) 1 



[J (b) + J 2 (b)} 



1 

A 



(39) 



b 2' 

where the first equality follows from the recurrence relation between Bessel functions J n . 
The absolute minimum of the sum of Bessel functions [Jo(b) + ^(fr)] equals —0.1323. The 
necessary condition for the disappearance of damping has thus the form 
2 



.4 



< 0.1323, or A = li^no [1 + H[l)] > A* = 15.1171... — 15.12 



(40) 



The above inequality shows that the nature of correlations could change for sufficiently high 
values of the surface fraction £ [Eq. fl29]) ] occupied by hard disks. The determination of the 
value of £* corresponding to the equality 

A* = 8C [1 + H(x = 1; O] = 15 - 12 

requires the knowledge of the contact value H(l) as function of the surface fraction. We 
have studied this question numerically. Accurate estimation of the precise value of £* cor- 
responding to A* is difficult because as we approach £* the iteration procedure demands an 
increasingly larger number of iterations to converge to a solution with the required accu- 
racy. Additionally, a computational domain over which the solution is sought must also be 
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progressively extended as we approach £*, since the decay of H(x) is very weak there. We 
estimated the value of by calculating A(£) for several values of £ in the range 0.5 < £ < 0.6 
and then extrapolating to larger values of £, as illustrated in Fig. [5j In this way, we obtain 
the estimate of £* m 0.622. 

Let us close this section with a comment on the ring approximation. In order to find the 
asymptotic behavior of correlations satisfying equation (1281) it is sufficient to consider the 
linearized version of equation ( 135]) . But this is equivalent to the replacement of the factor 
A = 2nna 2 [1 + H(l)] by 2nna 2 = 8£. The inequality (jlQJ) becomes then 

8£ > 15.12, or £ > 1.89 (41) 

representing for £ a physically impossible condition (beyond close packing). The ring ap- 
proximation is unable to describe a qualitative change in the hard disk correlation function. 
For any accessible surface fraction it predicts exponentially damped oscillations. 

C. Structural transition in a hard sphere fluid 



The integral equation satisfied by the correlation function H(x) of a three-dimensional 
hard sphere fluid for x 3> 1 has been analyzed in Ref. 15 within the ring approximation. 
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FIG. 5. The extrapolation of -A(£) dependence. The points correspond to the values of A obtained 
from the iterative procedure, the dashed line is given by A = A* ~ 15.12 

The frequencies of the exponential modes exp(/tx) describing the long distance decay of 
correlations are in this case solutions of the equation (see eq. (22) in^) 

k 2 = Anrra 3 ( ^ - ch^ (42) 



K 

In passing to the Kirkwood superposition approximation we need only to replace the factor 
Anna 3 in the above equation by A = 4nira 3 [l + H(l)]. This fact follows from the previously 
made remark that the ring approximation represents exactly the linearized Kirkwood theory. 
In order to determine the range of values of A > beyond which the exponential damping 
becomes impossible we explore the possibility of vanishing of the real part a of the complex 
number k = a + ib. Equation ( 1421) reduces then to 

,9 , f , sin6\ , b 3 

6 = A cos 6 — , or A = 43 

\ b J b cos b — sin b 

The absolute minimum of the function 

V(b) = 7 ^—— h (44) 

b cos b — sin b 

in the region where y(b) > equals 

y min = 34.81... (45) 
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Hence, for A > 34.81 equation (142!) acquires purely imaginary solutions and the exponential 
damping vanishes. It is quite remarkable that sixty years ago J.G. Kirkwood, E.K. Maum 
and B.J. Alder- concluded from numerical studies of the integral equation for H(x) that 
when A exeeds 34.8 the correlation function H(x) is not integrable any more. Our method 
provides a simple analytic confirmation of this result. 



D. The hard rod fluid: testing the method 

The rigorous calculation of the two-particle correlation function for hard rods (see e.g. 
Ref. 17) shows that its structure corresponds to exponentially damped oscillations at all 
possible densities. One finds 

na[H{x) + I] = (jr9{x-(k + 1)] ^ - ^ + exp{-C[a - (k + 1)]} (46) 

where 

na 



c 



1 — na 

In fact, the superposition law turns out to be exact for a one-dimensional hard rod fluid^, 
and the second equation of the equilibrium hierarchy takes a particularly simple form 

H\x) = na[6(x - 2)H(x - 1) - H(x) - 6(2 - x)] [H(l) + 1] (47) 

The formula ( 146]) represents the solution of (|4"T|) . When x > 2, we find 

H'(x) = na[H(x - 1) - H(x)\ [H(l) + 1] (48) 

Applying the same method as that used for hard disks we look for exponential modes exp(nx) 
solving ( 14H1) . The complex frequency k satisfies the equation 

« = na[H{l) + l]{e- K - 1) = C(e~ K - 1) (49) 

It turns out that all solutions of equation (T49l) can be expressed in terms of the multivalued 
Lambert W function. Indeed, the special function W(z) is defined on the complex plane by 
the equation 

z = W(z) exp[W(z)\ (50) 

But (jUJ) can be rewritten as 

(« + C) exp(/c + C) = CexpC (51) 
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FIG. 6. The real part of n m = —Q + W m (Cexp£), m=l,. . . ,5 (top to bottom) as a function of 



The principal branch of Lambert function, Wq(z), obeys Wo(xe x ) = x for real x, thus k = 
for that branch. However, other branches, W m (z) with m > 0, give values of K m with a 
negative real part, corresponding to the decay of the the correlation function. The first five 
solutions, K m ((), m = 1, . . . , 5, are presented in Fig. [6] as functions of no = C/0- + C) (modes 
with larger m decay faster). 

The negative real part of K m , m > is different from zero for any accessible density, and 
vanishes only at close packing where na — 1. There are thus no purely imaginary solutions 
k = ib of equation f )49|) . This confirms the correctness of the method, showing the absence 
of any structural transition in one dimension over the whole range of densities below close 
packing. 

As illustrated in Figs. [7] and [HI the correlation function H(x) rapidly approaches the 
asymptotic form given by the slowest decaying mode, Ae KlX . Note that not only the ex- 
ponential decay rate but also the oscillation period of the function agree with that given 
by Ae KlX starting from x ~ 3a. This shows that the other modes play a negligible role in 
influencing the behavior of H(x) for intermediate and large x values, thus lending further 



na = C/(l + C)- 



It follows that 



K = -( + W((exp() 



(52) 
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support to our approach of focusing on the slowest decay mode only. 

VI. DISCUSSION AND CONCLUSIONS 

The question of whether a system of particles interacting via a purely repulsive potential 
(only) can undergo a phase transition has been under continuous investigation since first 
posed and addressed by Kirkwood over 70 years ago. For a system of hard disks, the 
first numerical evidence was provided by Alder and coworkers^ 1 ^. The data reported in 
Ref. 21 showed that the hard disk freezing transition occurred at a density smaller than 
the density of closest packing [corresponding to an area fraction of £ = 7r/\/l2 = 0.90690], 
and suggested that the liquid to solid transition was first order. The most recent Monte 
Carlo simulations (on a system of 4 ■ 10 6 disks) of Mak 16 suggest that melting consists of 
a continuous transition from the ordered solid to an intermediate (hexatic) phased"— at a 
packing fraction rj = 0.723, and either a very weak first-order or a continuous transition 
from the intermediate phase to the fluid phase at a packing fraction rj = 0.699. 

From the analysis and numerical evidence presented in Section 4.B, we have calculated 
the area fraction £ at which a structural change in the hard disk fluid can take place, viz., 

17 



H 



0.1 



1 




A 




0.01 



0.001 



io- 4 

10~ 5 



x 



2 



4 



6 



8 



10 



12 



14 



FIG. 8. The absolute value of H[x) — 1 for hard rod fluid (solid) at no = 0.8 and the exponential 
asymptote H(x) = Ae Re ^ x 

~ 0.622. Converting this area fraction to a packing fraction gives rj* = £*/£o ~ 0.718. 
When compared to the estimate reported by Mak for the transition from the ordered solid 
phase to a dense fluid phase, r\ ~ 0.723, one finds the two values are in substantial agreement. 
Also of interest is our prediction of a structural transition in a system of hard spheres. As 
noted in Section 4.C, Kirkwood, Maum and Alder- found that for values of A > 34.8, no 
solutions of the YBG and Kirkwood integral equations exist for which, in their notation, 
x 2 [g (x) — 1] is integrable. We find that the value of A beyond which exponential damping 
becomes impossible is A = 34.81. Hence, the results for hard spheres appear to be in exact 
agreement. When we consider a system of hard rods, the analytic method developed here 
confirms the absence of any structural transition in d=l for the entire range of densities 
below close packing (a result which was already known to Rayleigh^ and Korteweg^). 

Our prediction of a structural phase transition is based on the analysis of an integral 
equation whose derivation assumed sufficiently fast decay of correlations. If a phase becomes 
ordered and correlations do not decay, the integral equation to which our method was applied 
does not hold. This can preclude the possibility of studying a region where a new equilibrium 
phase may be formed. In particular, the identification of a region intermediate between 
melting and freezing (see following paragraph) and the characterization of an ordered (solid) 



18 



phase is certainly beyond the scope of our approach. To study these questions within the 
Kirkwood superposition approximation, one must go back to the original YBG hierarchy 
equation (Kirkwood's closure does not assume the rapid decay of correlations). 

To develop this point further, there are three structural aspects of the hard disk transition 
that are not captured by the method developed in this paper. First, we find no evidence for 
the existence of an intermediate, or hexatic, phase predicted by the Kosterlitz, Thouless, 
Halperin, Nelson and Young (KTHNY) theory of d=2 melting^—, and supported by Mak's 
Monte Carlo simulations. Second, we find no evidence for the development of a shoulder on 
the second maximum of disk radial distribution function in the vicinity of the freezing tran- 
sition (rj = 0.686), reported by Truskett et al.— based on molecular dynamics simulations, 
and later correlated with structural rearrangements occurring at increasing disk density^ - —. 
Thirdly, we cannot confirm the existence of regions of five-fold coordination in the dense 
fluid phase, first predicted by Bernal^"— based on his ball and spoke model of a random 
assembly of hard-core particles, although it has been conjectured that the hexatic phase 
might be correlated with randomly dispersed regions of 5-, 6-, and 7- member disk clusters 
forming percolated tessellations that span the d=2 spaced. 

The larger point, however, relates to the original Kirkwood prediction, viz. that a system 
of particles interacting via purely repulsive forces, here hard disks but also hard spheres, 
can undergo a phase transition. Although not widely accepted at first, following the work 
of Onsager on the isotropic-nematic transition in a d=3 dimensional system of thin hard 
rods^, there developed a gradual realization that a phase transition can be entropy driven. 
As elaborated by Frenkel^, in hard-core systems the entropy in the ordered phase is actually 
larger than the fluid phase; quoting directly, "the entropy decreases because the density is no 
longer uniform in orientation or position, but the entropy increases because the free-volume 
per particle is larger in the ordered than in the disordered phase." The present contribution 
provides further evidence for the essential correctness of Kirkwood's insight. 

Appendix A: Derivation of Eq.(5) from the BBGKY hierarchy 

Consider a gas of hard disks of mass m and diameter a. We denote by j = (r^v,), 
j = 1,2, ... the one-particle state in which a disk has position r^ and velocity Vj. The average 
number density of s-particle clusters occupying at time t the s-particle state (1,2, ...,s) is 
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called the s-particle reduced distribution / s (l, 2, s; t). 

The dynamical evolution of the hard disk fluid is described in the thermodynamic limit 
by the BBGKY hierarchy equations. The second of them establishes a relation between / 2 
and / 3 

(§t + vi ■ ik. + v2 ■ ik ~ 2) ) /a(i ' 2]t)= I d3 3) +T(2 ' 3))/3(i ' 2 ' 3; t] (A1) 

The effects of binary collisions are described by the operator T(i,j) 

T(i,j) = a J dfrv 12 • <r0(v y ■ &) [<J(r - - <r)b a - <J(r y + a)\ , (A2) 

Here the Dirac (^-distributions restrict the distances |rjj| between the centers of the disks 
at the moment of impact to their diameter a = \cr\. The vector cr = a a is oriented 
perpendicularly to the surface of colliding disks at the point of impact. The action of the 
operator b a consists in replacing the velocities Vj,v, by their precollisional values v^,v^ 
corresponding to the inverse elastic collision. As in elastic collisions the kinetic energy is 
conserved, in the case of products of Maxwell distributions 

( 777- i 

we find 

b&MvMvj)] = MMvj) = <Kvi)<Kvj) (A3) 
Hence, in the case of equilibrium reduced distributions 

/ s (l, 2, s) = n s ( ri , r 2 ...r s )0( Vl )0(v 2 )...0(v s ) (A4) 
equation (lAip takes the form 

( Vl ' lL~ + V2 ' lL ° J daVl2 ' a5 ( Vi d ~ ra 2(ri 2 )0(vi)0(v 2 ) = (A5) 

a J dr 3 dv 3 / da [vi 3 ■ a 5{r u - a) + v 23 ■ a <5(r 23 - cr)] 

xn 3 (r 12 , ri 3 , r 23 )0(vi)0(v 2 )0(v 3 ) 

We can divide both sides of equation (1A5j) by the product 0(vi)0(v 2 ), and perform integra- 
tion over the v 3 variable. Moreover, we introduce explicitly the excluded volume factors by 
using equations (JT]) , (SJ) - The hierarchy equation becomes 



^v i2 • t^- - a J dav 12 - a5(r 12 - crU 0{r X2 - o)y 2 (r l2 ) 



(A6) 
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na9{r l2 - a) j dr 3 J da [vi • & 8(r 13 - a)9(r 23 - a) - a)y 3 (r 12 , cr, r 23 ) 
+v 2 • <x<5(r 23 - a)9(r 13 - a)y 3 {r 12 , r 13 , a)] 

We now use the identity 

d f 
V12 • ^ — 9(r 12 - cr) = cr / d&vi2- &5(v 12 - a) (A7) 
Or 12 J 

which reduces the left hand side of (IA6I) to 

d 

L = 9{r X 2 - cr)vi2 • ^ — 2/2(^12) (A8) 
Or\2 

On the right-hand side owing to the presence of ^-distributions we can perform integration 
over variable r 3 thus obtaining 

R = na9 (7-12 - cr) y d(TVi2 • <x0(|ri 2 - a) \ - a)y 3 {r X2 , a, |r 12 - cr) |) (A9) 

As the equality L = R must hold for any value of the relative velocity vi 2 , we finally find 
(when ri2 > cr) 

dV ^ 12 =na J dar 12 ■ a9(\r 12 - a\ - a)y 3 (r 12 ,a,\r 12 - a\) (A10) 
which is equation ([5]) of section 2. 

Appendix B: Derivation of Eq.(27) 

We perform here the angular integration in the contribution to the correlation function 

/*00 />27T 

I**(x) = —— / dz dd> cos <j>6{z - 2 cos 4>) H (v^ 2 + 1 - 2,2 cos 0) (Bl) 
2 Jx Jo 

poo P n /2 

— — I dz d(p cos (p[9(z — 2 cos 4>)H(^/z 2 + 1 - 2z cos 0) - i/(Vz 2 + l + 2zcos0)] 

Jx Jo 

Changing the order of integrations with the use of the asymptotic decay of the correlation 
function we arrive at a convenient formula 

/■tt/2 rx+cos<j> j 

I**(x) = - d0cos<W ds 9 <(s- cos 0) H '(yj s 2 + sin 2 0) (B2) 

JO J X — COS0 

Putting then /i = sin we get 

= - j dfi j ds9{x + - s)9{s - x + v 7 ! - /W 2 ) (B3) 
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We now introduce a new integration variable 

z = \J s 2 + /i 2 

As zdz = sds we find 

I**(x) = - [ dn [ dz Z =H(z)9(^z 2 - fi 2 -^l^)9( y ^Jl 2 -\x- v / z 2 - fi 2 \) 

Jo J V Z 2 — fi 2 

= - I dz zH(z)6(z - 1) I d/i 1 9(y/l -ii 2 - \x- \/z 2 -n 2 \) 
J Jo v z 2 — /i 2 

Here the inequality y/l — \i 2 > \x — \J z 2 — /i 2 | is equivalent to 

/-= ■= x 2 + z 2 -1 

— M > ^ 

which leads to the formula 

,.2 i ^2 



r*(x) = - / cfe ztf - 1) / dp ^—J ){^^ 2 - X + * x * ) (B5) 







used in equation (|Tj 
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(B4) 



Putting fi = zu we find 

= - dz zH(z)6(z - 1) / A/ ; gfy/T^ 2 - — ) (B6) 

J Jo VI- v 2 %xz 

The change of the integration variable w = \J\ — v 2 yields the formula 

r r 1 dw r 2 + z 2 — i 

I~(x) = ~ dz zH{z)6{z - 1) / "7=^ " ^ ) ( B? ) 

= - / ^ - 1) r -^=s«a - x2+ /~ 1 ) 

= ~ J dz zH(z)6(z - 1)0[1 - (x - ^) 2 ] 1 1 - arcsin ^ 
Finally, we arrive at the result 

r x+1 (x 2 + z 2 -\\ 

I**( x ) = - dz zY{z)6{z - 1) arccos —- (B8) 

Jx-i V %xz / 
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